Stein, Belgium is an agriculture village which have many problems with soil pollution. Due to the situation, the authorities created the Meuse dataset that gives information about the top soil heavy metal concentrations (ppm), along with a number of soil and landscape variables. Heavy metal concentrations are bulk sampled from an area of approximately 15 m x 15 m (Pebesma, 2020).

This dataset contains the following columns:

For interpolate zinc values we use ordinary kriging method. This is a gaussian approach provides always the best linear unbiased prediction for all predictions based on covariance/variogram estimates.

To use kriging two hypothesis must be fulfilled:

1. The first moment of the sample is stationary.

2. The correlation between random variables solely depends on the spatial distance and is independent of their location.

After estimate the mean on different spatial bootstrap samples (it is not shown here), we notice that all random variables follow a uniform distribution (first hypothesis). The no trend in the spatial autocorrelation (second hypothesis) of zinc values was measure using an empirical/experimental variogram that was fitted to a spherical variogram model (see left image). We notice that the range is stabilized at around 850 meters.

Using the spherical variogram model, the kriging weights were calculated by deriving minimum variance and substituting the semivariogram function into the kriging equation. The left map shows us the kriging prediction results and the local variance error for zinc values.

out-of-sample testing or cross-validation is a statistical method that could be used to estimate the local and global error of spatial models.

For this example we use 5-fold out-of-sample testing (see left image). In this algorithm, The original sample is randomly partitioned into k equal sized subsamples. Of the k subsamples, a single subsample is retained as the validation data for testing the model, and the remaining k − 1 subsamples are used as training data. The cross-validation process is then repeated k times, with each of the k subsamples used exactly once as the validation data. The k results can then be averaged to produce a single estimation.

The advantage of this method over repeated random sub-sampling (see right image) is that all observations are used for both training and validation and each observation is used for validation exactly once. The 5-fold cross-validation is used for zinc spatial prediction.

It can be seen that areas near the river present positive large residuals, while values in the inner area, at least most of it, present negative and close to zero values.

These results show that in future fieldworks, it would be ideal to collect more data from areas near the river rather than in inner areas since it would generate a lower global mean error in the spatial prediction.

If we random sampling with no replacement the meuse dataset and also considering different samples groups (25, 50, 75, 100, 122, and 155). We could indirectly measure the global spatial variance of our experiment. The left image shows the global interpolation error in the y-axis (RMSE) and different samples in x-axis.

To improve the stability and accuracy of our results, we use bootstrap aggregating (running 100 times the same experiment for each random sampling). The results reveal that the less number of points, the higher RMSE and lower stability.

One of the most important points to create a good spatial prediction model is determine the spatial dependence rule. In Inverse Distance Weighting (IDW), the power value is the most important parameter which have to be measured (e.g. with cross-validation). On the other hand, in ordinary kriging, we need to adjust the variogram and to be more specific the range, which is the maximum distance up to which the spatial autocorrelation exists.

In this experiment we run kriging ordinary, with different range values (200, 400, 600, 800, 1000 and 1200). The results illustrates that low range values generate a bull eye effect on the results while larger range values generate smoothed regions.

After analyzing the variogram, it can be noticed that zinc values present a clear spatial autocorrelation. In addition, the 5-fold cross-validation shows us that the global error of the spatial model proposed here (\(RMSE_{global}\) ) was 0.26 and this value could exponentially increase if the number of points is reduced.

These results can help decision-makers to design future fieldworks taking into account the error map (prior probability) generated in this report.